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Critical velocities for deflagration and detonation triggered by voids in a REBO high 

explosive 

S. Davis Herring, 1, 2 [j Timothy C. Germann/'LJ] and Niels Gr0nbech- Jensen 2 ||| 
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The effects of circular voids on the shock sensitivity of a two-dimensional model high explosive 
crystal are considered. We simulate a piston impact using molecular dynamics simulations with 
a Reactive Empirical Bond Order (REBO) model potential for a sub-micron, sub-ns exothermic 
reaction in a diatomic molecular solid. The probability of initiating chemical reactions is found to rise 
more suddenly with increasing piston velocity for larger voids that collapse more deterministically. 
A void with radius as small as 10 nm reduces the minimum initiating velocity by a factor of 4. 
The transition at larger velocities to detonation is studied in a micron-long sample with a single 
void (and its periodic images). The reaction yield during the shock traversal increases rapidly with 
velocity, then becomes a prompt, reliable detonation. A void of radius 2.5 nm reduces the critical 
velocity by 10% from the perfect crystal. A Pop plot of the time-to-detonation at higher velocities 
shows a characteristic pressure dependence. 

PACS numbers: 62.50.Ef, 47.40.Nm, 82.40.Fp 
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Microscopic defects in solid high explosives can have 
dramatic effects on the sensitivity of the bulk material to 
initiation by shock waves [l|. When a shock that is in- 
sufficiently strong to ignite the material directly encoun- 
ters a defect, reflection and refraction redirect its energy. 
Where the energy density is reduced, little changes; the 
shock was already inert. However, a local, temporary in- 
crease in energy density may drive exothermic chemical 
reactions; the resulting hotspot will not generally cause 
detonation directly, but will increase the pressure behind 
the shock and thus its strength. As the shock increases 
in strength, more energy is available for focusing and less 
focusing is needed to initiate further reactions. The pos- 
itive feedback that results is the principal driver of the 
transition from shock to detonation in heterogeneous ex- 
plosives [2[. The details of the defect feedback process 
remain poorly understood [2| , and so practical questions 
like "To what extent would a population of voids (with 
some distributions in space and size) reduce the critical 
velocity of this explosive?" go unanswered. 

Spherical voids are a common defect in cast and formed 
explosives Q . When a sufficiently strong shock wave en- 
counters a void, the leading surface is ejected into the 
void, and the resulting gas is compressed as the void 
collapses; jetting may also occur [l], 0, [4J. The shock's 
energy is focused onto the downstream pole of the void. 
If the void is large enough and the shock strong enough, 
chemical reactions result. This process is illustrated in 
Fig.ffl 
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FIG. 1: (Color online) Snapshots from a 124 x 135 
atom low- velocity simulation with r = 20 A and u p 
Undisturbed material at left is the piston; green, red, and pur- 
ple atoms are products, radicals, and clusters respectively. 
(a) Just after impact, (b) 200 fs before the void finishes col- 
lapsing, (c) End of the simulation; shock has reached the free 
surface 1.17 ps after impact. 



Several attempts to identify the atomistic mechanism 
of chemical initiation in collapsing voids have been made 
using molecular dynamics (MD) simulations |5[. Mint- 
mire et al. [6| determined that energy was efficiently 
transferred into the molecular vibrational modes of a 
nonreactive diatomic molecular solid only when the col- 
lapse of the void was turbulent and involved the disinte- 
gration of its walls; such vibrational excitation is thought 
to be a necessary precursor to chemical reaction. White 
et al. 7] found that randomly placed circular voids signifi- 
cantly affected the response of ozone crystals under shock 
loading. Germann et al. [8[ observed that reactions oc- 
curred some time after the ejecta collided with the down- 
stream wall (and, with a periodic array of voids, could 
lead to detonation), and that the reduction in critical 
shock strength for ignition depended on the orientation 
as well as the size of elliptical voids. In particular, sensi- 
tivity was observed to increase with the width of a pla- 
nar gap, suggesting that heating via recompression of the 



ejected gas was important for initiation. Holian et al. [9| 
extended the planar gap analysis with a Lennard- Jones 
potential and derived an expression for the temperature 
increase due to recompression. 

Hatano [10| considered the non-equilibrium mechan- 
ics of the ejected material in cuboidal voids and ob- 
served that the frequency of energetic molecular collisions 
reached a maximum after temperature did and could vary 
independently of that maximum. Shi and Brenner jllfl 
considered infinite rectangular voids in Ns cubane crys- 
tals and observed reactions following almost immediately 
after the initial downstream jet impact, with individual 
molecular collisions at the impact point leading directly 
to dissociation. Turbulent destruction of the void walls 
and focusing of the ejecta by the walls were observed 
to increase sensitivity in that system, but recompression 
after the jet impact (or in the absence of any jet) was 
not. 

To better understand the output of void-based 
hotspots (as a function of void size and input shock 
strength) and their ability to precipitate detonation, we 
use MD simulations of supported shocks in an assortment 
of two-dimensional samples having one circular void (or, 
more precisely, a periodic line of voids) each. In this 
initial exploration we consider only one generation of 
hotspots, excluding the feedback process between a shock 
and the sequence of voids it encounters. For high enough 
shock strengths, however, a similar feedback can develop 
among scattered reactions in the shocked material. A 
similar transition to detonation thus occurs nonetheless, 
as depicted in Fig. [5J 



II. METHOD 
A. Approach 

We perform two separate investigations in different ve- 
locity regimes to study the transitions, as radius and/or 
velocity increase, from no reactions at all to slow defla- 
gration and then to a detonation wave. 

Low-strength shocks may collapse a void without pro- 
ducing any chemical reactions. We simulate low-velocity 
pistons impacting samples with voids of various radii; 
Fig. [1] is taken from one such simulation. In each case, 
several initial conditions with varied random thermal ve- 
locities are considered so as to obtain a probability of 
initiation. Very small voids (and the special case of zero 
void size) are excluded from this study because sponta- 
neous reactions would compete with those triggered by 
the void and spoil the results. 

At higher impact velocities, the hotspot always reacts 
but may produce a detonation wave only much later if at 
all. It is difficult to establish with a finite sample that a 
detonation would never develop from an observed defla- 
gration, but the progress of the deflagration at moderate 
velocities and the promptness of the detonation at high 
velocities may be used to bracket the critical velocity. 



The high-velocity piston impact simulations are similar 
to those in the low-velocity study: they involve various 
void radii and start from several thermally random initial 
conditions for each case. Figure [5] is taken from one such 
simulation in which the sample detonated. 

The stochastic initiation process is studied in samples 
having voids of radius 1 to 10 nm with piston veloci- 
ties of up to 3.5 km/s. The high-velocity study of the 
transition to detonation involves samples having voids 
of radius 1, 2.5, or 5 nm as well as the perfect-crystal 
case (r — 0). There, the pistons have velocities of 
Up = 2.95-4.90 km/s and produce pressures of 10.8- 
20.9 N/m. Such two-dimensional pressures are difficult 
to interpret physically, but following [12j we may suppose 
that 1 N/m sw 2.5 GPa. 



B. Model 

The Reactive Empirical Bond Order (REBO) "AB" 
potential (originally developed in [5|, [l2|, [l3|]) describes 
an exothermic 2AB -> A2 + B2 reaction in a diatomic 
molecular solid and exhibits typical detonation proper- 
ties but with a sub-micron, sub-ns reaction zone that 
is amenable to MD space and time scales. Heim et 
al. modified it to give a more molecular (and less plas- 
malikc) Chapman- Jouguet state [lj]- We utilize the 
SPaSM (Scalable Parallel Short-range Molecular dynam- 
ics) code [15] and that modified REBO potential ("Mod- 
ellV"). A standard leapfrog- Ver let integrator is used with 
a fixed timestep of 0.509 fs in the NVE ensemble. 

Each two-dimensional sample is a rectangle of herring- 
bone crystal with two AB molecules in each 6.19 x 4.21 A 
unit cell. The shock propagates in the +z direction; the 
samples are periodic in the transverse x direction. A cir- 
cular void is created by removing all dimers whose mid- 
points lie within a circle of a given radius (see Fig. [Ha)). 
In the low- velocity cases, the sample is made large enough 
to prevent interaction between the periodic images of the 
void until the collapse is finished and the material has or 
has not reacted. Depending on the size of the void, 672- 
49700 atoms are simulated. The high-velocity samples 
are 1 um (2381 unit cells) long and two void diameters 
wide (or 10 nm in the case of no void), with the void 
center four diameters from the piston face; they include 
66622-313086 atoms. 

Three layers of unit cells at the — z end are frozen 
to serve as a piston (of infinite mass), and the rest 
is assigned a tiny but finite temperature (low veloc- 
ity: 11.6 mK, high velocity: 1.00 mK) and a bulk velocity 
v z = —Up directed into the piston. (The temperatures 
must be small to avoid melting the material: the 5 meV 
depth of the van der Waals well corresponds to just 58 K.) 
To reduce the correlation between the different histories, 
the initial thermal velocities are allowed to thermalize 
for 1 ps (5 ps for the low-velocity study) before the bulk 
velocity is applied. 



FIG. 2: (Color online) Snapshots of the forward ~40% of a 204 x 10013 A, 313086-atom high-velocity simulation with r = 50 A 
and Up ~ 4 km/s; colors as in Fig. \T\ (a) The periodic images of the hotspot have merged, (b) Reactions develop between 
the deflagration and the shock front, (c) Detonation has commenced, separated from the original deflagration; the zone of 
increased dissociation near the shock corresponds to the overdriving in Fig. [7] 



C. Analysis 

We define two atoms as bound if one does not have 
escape velocity with respect to the other (taking the po- 
tential's repulsive core into account 14]). At regular in- 
tervals, the atoms to which each atom is bound are noted. 
Two atoms each bound only to the other are deemed a 
molecule; heteronuclear AB molecules are reactants and 
homonuclear molecules are products. Many of the results 
derive from a count of such reacted molecules. Atoms 
bound but not in a molecule are termed clusters; Fig.QJc) 
involves all four possible labels. 

To identify in which simulations and at which times 
detonations develop, we measure the position of the 
shock wave (whether reacting or not) at intervals of 
51 fs throughout each simulation. The shock positions 
are obtained by finding the pair of adjacent columns of 
well-populated computational cells (of thickness Az ps 
0.53 nm) with the largest difference in (v z ). The identi- 
fied positions are thus only accurate to Az and are occa- 
sionally much too small (when some local fluctuation in 
the shocked region is misidentified as the shock). Before 
seeking the detonation transition, the positions are fil- 
tered by removing all values smaller than any preceding 
them. 

Detonation transition times are extracted from the fil- 
tered shock positions by finding (Obi, z i), (^2, £2), Ots, 23)) 
triples with Z3 — Z2 and zi — zi each > 10 nm that max- 
imize the weighted second derivative 
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where the additional ^/m2 slope factor favors the tran- 
sition to detonation over any sudden acceleration asso- 
ciated with mere deflagration. Maximum scores greater 
than a manually chosen threshold of 3.9 x 10 12 m 3 / 2 s~ 5 / 2 
are taken to indicate detonation transitions. (At low ve- 
locities, the void collapse can be mistaken for a transition. 
These false positives are easily identified by the large gap 
between them and the true detonations.) 

We also determine the shock pressure p associated with 
a piston velocity u p from the shock position data. We 
measure the shock velocities u s from the beginnings of 
voidless simulations (to minimize the effects of reactions 




FIG. 3: Hugoniot for the unreacted material for piston veloc- 
ities through the largest used in the high-velocity study. The 
line is a smoothing quadratic fit through the data points in 
its interval. 



on the shock positions) and calculate p — pqu s u p (the 
Hugoniot jump condition with po — 0). For reference, 
the resulting unreacted Hugoniot is shown in Fig. [J2 the 
pressures used in Section flll Bl are smoothed by sampling 
from a quadratic fit to the highest-velocity data. (The 
product Hugoniot is presented in [14j.) 



III. RESULTS 

A. Low-velocity regime: probability of initiation 

Initiation probabilities were derived from at least 20 
realizations of each of 1560 radius-velocity pairs (73553 
simulations total). The probabilities obtained for three 
radii (the largest, the smallest, and an intermediate value 
with high-quality data) are given in Fig. 0J At each ra- 
dius (10, 12, 14, . . . , 98 A), the critical velocity u p ^ 50 
(as explored in Q) and transition sharpness a were de- 
termined by scaling and shifting the sigmoid function 
Pq(x) := (1 + erf x)/2 in velocity to fit the measured 





FIG. 4: Probabilities of initiation as a function of piston ve- 
locity at three selected void radii, each with a sigmoid fit. 
Larger voids shift the transition to lower velocities and make 
it sharper. At r = 20 A, 1000 simulations were performed per 
velocity, so the statistical noise is much reduced. 



FIG. 5: Fit parameters [see Eq. ©] for the initiation prob- 
ability as a function of piston velocity at each void radius, 
with fits to them described in the text. Larger inverse widths 
correspond to more sudden transitions from no reactions to 
guaranteed reactions. 



probabilities: 



P(u p ) 



1 + crf a(u p - u p , 50 ) 



(2) 



(Other similar functions [e.g., the logistic function 
L(x) := (l + e _x ) -1 ] were considered; Po(x) was selected 
because its width-slope product was judged most similar 
to that in the data.) Figure H] also contains the sigmoid 
fits for its three radii. 

The parameters in Eq. ((2|) for all radii are given in 
Fig. [5] Also shown are fits to the center (the 50% con- 
tour) and scale parameters as functions of radius: 
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and 
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where a = 1.99 x 1CT 6 m 2 /s, (3 = 718 m/s is the asymp- 
totic value for large voids, 7 = 3.88 x 1CT 4 s/m, and 
6 = 0.628. The critical velocity for initiation decreases 
with radius, as the void focuses more of the shock and cre- 
ates higher temperatures @,Q; at r = 10 nm it is a factor 
of 4 lower than the minimum velocity observed to initi- 
ate reactions in a voidless sample. (Several other models 
were considered for the critical velocity, but Eq. ^ was 
far more successful than the other fits.) The width (in 
velocity) of the transition also decreases with increas- 
ing void radius, as the hotspot development becomes 
less dependent on the stochastic behavior of individual 
molecules. 



B. High-velocity regime: transition to detonation 

Each of 79 radius-velocity pairs was simulated 20 
times, each until the shock broke out at the free sur- 
face. In each simulation, the maximum number of re- 
acted atoms observed (typically also the final count) was 
noted. The average count for each r-u p pair, as a frac- 
tion of the total number of atoms, is given in Fig. [6] Even 
when the reaction consumes the entire sample, the con- 
version to A2 and B2 is not complete; reaction extents 
> 70% indicate reaction of the entire sample (and thus 
detonation, since the reaction reached the shock front). 
At velocities above 4.4 km/s, detonation is inevitable re- 
gardless of void size (or even presence) ; the reaction ex- 
tent decreases slightly as velocity increases further be- 
cause the equilibrium at the higher energy states is less 
completely reacted. 

At velocities of 3.0-3.6 km/s, the hotspot consistently 
establishes a growing deflagration in the sample that 
merges with its periodic images and becomes planar but 
is left behind by the shock and does not become a det- 
onation. This process produces the consistent reaction 
extents of approximately 20% that appear for all non- 
zero void sizes in Fig. [B] The drop-off in the r — 10 A 
data below u p = 3.4 km/s is largely due to failure to 
create a reacting hotspot, rather than due to reacting 
hotspots quenching. For larger voids, that drop-off moves 
off below the u p range of the plot, and the "shoulder" of 
deflagration widens. The step up to detonation moves 
more slowly and is of limited utility in identifying a crit- 
ical velocity because more simulation time (with a larger 
sample) might allow some of the deflagrations to become 
detonations. 

The shock positions from a representative simulation 
that developed a detonation are shown in Fig. [7] (the 
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FIG. 6: Average maximum extent of reaction at each pis- 
ton velocity and void radius. Values in the plateau indicate 
detonation. 
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FIG. 7: History of the shock position in Fig. [2]s simulation. 
The times of those snapshots are marked. The detonation 
transition is clearly visible; the dotted straight lines serve to 
illustrate the shock's acceleration before and after it. 



Az-scale roughness is invisibly small). Prior to the tran- 
sition, reactions developing behind the shock accelerate 
it; this acceleration occurs even in the samples with- 
out voids, where one would expect the shock to accel- 
erate only at transition (when the homogeneous initi- 
ation catches up from the piston face). The conditions 
that inspire a detonation within such small homogeneous 
samples entail a chemical induction time so short that 
randomly distributed hotspots appear spontaneously and 
drive the detonation in the same fashion as in the het- 
erogeneous case. 

Temperature profiles at four times during the same 
simulation are given in Fig. [3] The shock positions, the 
spread of several deflagrations throughout the shocked 
material, and the transition to detonation at z — 290 nm 
are evident. The persistent dip and spike astride the 
transition point seem to indicate that, when the reaction 
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FIG. 8: Temperature profiles from Fig. [2js simulation over 
approximately the same region of the simulation. The first 
three traces correspond to those snapshots; the last is from 
the end of the simulation and illustrates the stability of the 
features. Temperatures of 2.5, f2.5, and 15 kK correspond 
to the inert shock, deflagration in the shocked region, and 
detonation respectively. 



zone was just behind the shock but had not yet merged 
with it, the energy from the reaction went to accelerating 
the shock (and thus further heating the material being 
shocked) rather than to heating. At the latest time, the 
boundary between the regions of deflagration and det- 
onation remains well-defined, although the temperature 
spike has advected and diffused somewhat. 

The minimum transition times for each r-u p pair (for 
which any detonations were observed) are given in Fig. [5] 
For comparison, the medians are also given for those 
voidless velocities for which a majority of simulations 
produced detonations; the medians for other void sizes 
behave similarly, with a nearly constant ratio between 
the minimum and median times. We see pressure de- 
pendences of p~ 13 - 77 with no void and p~ 9 - 95 with any 
size of void, both much larger than the (space-pressure) 
exponents of -1.6 for PBX-9501 [jj] and -4.5 for PBX- 
9502 [17] . This discrepancy may arise from the dimen- 
sionality of the system and the associated unusual pres- 
sure units. 

At very high pressures, the void is irrelevant even to 
the promptness of the detonation, as many reactions are 
initiated directly upon the piston face. (The largest voids 
even retard the high-pressure transitions, perhaps due to 
their greater distance from the piston.) At lower pres- 
sures, the presence of a void greatly accelerates the devel- 
opment of a detonation by providing a guaranteed source 
of significant deflagration, but the size of the void seems 
not to significantly affect the subsequent positive feed- 
back that develops the detonation. However, the tran- 
sition times for the three non-zero radii separate at low 
pressures, where the energy available from the hotspot 
becomes the determining factor in developing a detona- 
tion. It happens that this change in the pressure expo- 
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FIG. 9: (Color online) Pop plot of the minimum times to 
detonation, with power-law fits to the no-void data and to all 
the data with voids. The median times to detonation are also 
presented for the no- void case for illustration. 



nent (for each radius) occurs just as the transition times 
become longer than the simulation, and so corresponds 
precisely to the lower limit of detonations in Fig. [51 



IV. DISCUSSION 

We have obtained a simple form [Eq. ([3]), Fig. [5] for the 
piston velocity needed to create a reacting hotspot from 
a void of a given radius in a two-dimensional ModellV 
REBO high explosive crystal. In the limit of large voids 
it appears sufficient to give each dimer just 70 meV of 
kinetic energy (7% of the height of the repulsive core). 
A difference of less than 400 m/s in piston velocity is 
expected to switch from 10% to 90% chance of ignition 
for voids with r > 10 nm. 

At higher velocities, we have observed a transition from 
steady deflagration to steady detonation in the range 
of 4.0 ± 0.4 km/s (an average kinetic energy of 1.1 eV 
per atom). Its precise, radius-dependent location is non- 



trivial to establish, but these results suggest the form 
u c (r) m (3.8 + 0.34e~ r / 22A ) km/s. As larger samples 
(and thus longer simulation times) might allow more det- 
onations to finish developing, this expression is probably 
an overestimate. 

The detonation initiation mechanism observed is nei- 
ther the superdetonation associated with homogeneous 
explosives [2j nor purely the ignition and growth of 
(the single rank of) hotspots, although such growth is 
observed (see Fig. Ufa)). Rather, the heat and pres- 
sure produced by the deflagration outpace it, encourage 
further reactions throughout the shocked material (see 
Fig. H[b)), an d strengthen the shock until it ignites the 
material directly. (The sample width [or distance be- 
tween periodic images of the void] may affect the feed- 
back; a larger sample provides more possible reaction 
sites but also more material over which to disperse the 
void's output.) Extensions to this work may addition- 
ally consider the feedback resulting from the effects of 
a hotspot on further voids downstream. However, these 
results already suggest that even isolated nanoscopic fea- 
tures may directly affect bulk material behavior through 
the positive feedback inherent to energetic materials. 

We have also observed the typical power-law depen- 
dence of detonation induction time on (two-dimensional) 
shock pressure, although the exponent seems to change 
at lower pressures. To our knowledge, such a result has 
not previously been reported for an MD simulation. It 
remains to be seen whether such a relationship holds 
in three-dimensional systems (whose pressures are more 
meaningful) . 
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